/*
	ParticleSet: A class for representing the set of particles used in error correcting the Level Set
	Inputs: Size of grid and cell

	This class only stores a list with pointers to all the particles within the grid. Note that the
	particles are not stored per cell, but for the entire grid.

	Functions:
	Update		- takes as input a velocity grid and a timestep. Calls the update function for
				  each partricle and removes it if the particle has exited the grid.
	Resample	- Updates the radius for each particle. Only use this function is necessary
	Reseed		- Deletes all particles and creates new ones. Only use this function when
			      absolutely necessary

	Created by Emud Mokhberi: UCLA : 09/04/04
*/

#ifndef PARTICLESET_H
#define PARTICLESET_H

#include "main.h"
#include "vector.h"

#include "levelSet.h"
#include "particle.h"
#include "velocity.h"

class ParticleSet
{

    public:

      // Common to other classes - Parameters of the constructor

      int nx;
      int ny;
      int nz;
      double h;

      int size_vectors;

      double **gridUVW;
      double **gridCellUVW;
      int *flag;


      std::vector <int> *List_Border_Face_Air;
      std::vector <int> *List_Border_Face_Fluid;
      std::vector <int> *List_Border_Air;
      std::vector <int> *List_Border_Fluid;

      Scene *S;

      double* gridPhi1;
      double* gridPhi2;

      // Computed from the arameters of the constructor
      int dj, dk, hInv;
      double h2;

      //Parameters changed at every iteration

      double* gridPhi;
      double dt;

      double *expected_time;

       inline void GIJK(int index, int& i, int& j, int& k)
      {
          k = int(index/dk);
          j =(index%dk)/dj;
          i = index%dj;
      }

private:
	list<Particle*> particles;

public:
      ParticleSet(int nx_par,int ny_par, int nz_par, Double h_par,int size_vectors_par,double **gridUVW_par, double **gridCellUVW_par, int *flag_par,
      std::vector <int> *List_Border_Face_Air_par, std::vector <int> *List_Border_Face_Fluid_par,std::vector <int> *List_Border_Air_par, std::vector <int> *List_Border_Fluid_par,
      Scene *S_par, double *gridPhi1_par,  double *gridPhi2_par, double *expected_time_par)
    : nx(nx_par), ny(ny_par), nz(nz_par), h(h_par), size_vectors(size_vectors_par), gridUVW(gridUVW_par), gridCellUVW(gridCellUVW_par), flag(flag_par),
      List_Border_Face_Air(List_Border_Face_Air_par), List_Border_Face_Fluid(List_Border_Face_Fluid_par),List_Border_Air(List_Border_Air_par),List_Border_Fluid(List_Border_Fluid_par),
      S(S_par),gridPhi1(gridPhi1_par),gridPhi2(gridPhi2_par),hInv(1./h_par),expected_time(expected_time_par)
     {
        gridPhi=gridPhi1;
        dj = nx + 4;
        dk = (ny + 4)*dj;
        h2 = h*h;
      }

	~ParticleSet()
	{   for(Iterator it = particles.begin(); it != particles.end(); it++) delete *it;
        particles.clear(); }

	typedef list<Particle*>::iterator Iterator;
	typedef list<Particle*>::const_iterator cIterator;
	cIterator begin() const { return particles.begin(); }
	cIterator end()   const { return particles.end(); }

	void Update()
	{
	    dt = S->delt;

        static Particle* p;
        static Vector pos;
		for(Iterator it = particles.begin(); it != particles.end();)
		{
			p = *it;
			p->GetPosition(pos);
            p->Update(gridUVW , dt, hInv , size_vectors);
			p->GetPosition(pos);

			if( pos[0] < 0 || pos[0] > nx+1 ||
                pos[1] < 0 || pos[1] > ny+1 ||
                pos[2] < 0 || pos[2] > nz+1) { it = particles.erase(it); delete p; }
			else it++;
		}
	}

	void Resample(const LevelSet& levelSet) // updates particle radii;
	{
		static Vector pos;
        Particle *p;
		for(Iterator it = particles.begin(); it != particles.end();)
		{
			p = *it;
			p->GetPosition(pos);
			if(p->SetRadius(levelSet.SAMPLEPHI(pos), hInv)) it++;
			else { it = particles.erase(it); delete p; }
		}
	}
	void Reseed(const LevelSet& levelSet)	 // deletes particles and creates new ones
	{
        static Double phi,exp_time,ppn;
        static bool reseed, reseed2;
		for(Iterator it = particles.begin(); it != particles.end(); it++) delete *it;
		particles.clear();

		//ofstream reseedLog("reseedLog.txt" , reseedLog.app);
                int itend = (int) List_Border_Fluid->size();
                for(int itfor=0; itfor<itend;itfor++)
                {
                    int index=(*List_Border_Fluid)[itfor];
                    int i,j,k; GIJK(index,i,j,k);
                    exp_time =expected_time[index];
                    cout << "exp_time " << expected_time[index] << endl;
		    //FOR_LS
		    reseed = false;
		    reseed2 = false;
		    //int index = i + j*dj + k*dk;
		    //reseedLog << "cell (" << i << " , " << j << " , " << k << ")" << endl;
                    for(int dx=0; dx <= 1; dx++) {
                       for(int dy=0; dy <= dj; dy +=dj) {
                           for(int dz=0; dz <= dk; dz +=dk) {
                               phi = abs(gridPhi[index + (dx+dy+dz)]);    
                               /*if(phi < RESEED_THRESHOLD) reseed = true;
                               if(phi < h) reseed2 = true;*/
                               //reseedLog << "  neighbour (" << dx << " , " << dy / dj << " , " << dz / dk << ") exp_time " << exp_time;
			      
                               if((exp_time > 0) && (exp_time<EXP_TIME_LIMIT))
                               {
                                    //reseedLog << " passed the reseed test;";
                                    if(phi < RESEED_THRESHOLD) reseed = true;
                                    if(phi < h) reseed2 = true;
                                }
                                //reseedLog << endl;
                            }
                        }
                    }
                    if(reseed2) ppn = PARTICLES_PER_INTERFACE_NODE;
                    else        ppn = PARTICLES_PER_NODE;
			if(reseed) {
				for(int x=0; x < ppn; x++) {
					Vector pos(Double(i) + RandomFloat(),
							   Double(j) + RandomFloat(),
                               Double(k) + RandomFloat());
					particles.push_back(new Particle(pos, levelSet.SAMPLEPHI(pos), hInv , dj , dk));
				}
			}
                    }
                    // END_FOR_THREE
                    //reseedLog.close();
	}
        
        void draw(const LevelSet& levelSet){
            /*Draw a representation of the particles. The idea is to draw, per grid cell, only the first particle
            which is distant enough from the interface of the fluid, discarding the others that lie on this cell.
            */

            cout << "Particles draw.........................................................................." << endl;
            glTranslatef(-1.0f , -1.0f , -1.0f);
            glScalef(2.0f , 2.0f , 2.0f);
            
            float distThreshold = 0.5 * h;
            list<Particle*>::iterator partIter;
            typedef map<Vector , bool> partMap;
            partMap mapPos2Id;
            partMap::iterator mapIt;

            for(partIter = particles.begin() ; partIter != particles.end() ; partIter++){
                //Finding the cell where the particle lie.
                Vector pos;
                (*partIter)->GetPosition(pos);
                Vector cell(int(pos.x()) , int(pos.y()) , int(pos.z()));

                mapIt = mapPos2Id.find(cell);
                if(mapIt == mapPos2Id.end() && levelSet.SAMPLEPHI(pos) >= distThreshold){
                    mapPos2Id[cell] = true;
                    (*partIter)->draw(h);
                }
            }
        }
};

#endif
